// Copyright 2018 Ulf Adams
//
// The contents of this file may be used under the terms of the Apache License,
// Version 2.0.
//
//    (See accompanying file LICENSE-Apache or copy at
//     http://www.apache.org/licenses/LICENSE-2.0)
//
// Alternatively, the contents of this file may be used under the terms of
// the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE-Boost or copy at
//     https://www.boost.org/LICENSE_1_0.txt)
//
// Unless required by applicable law or agreed to in writing, this software
// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied.

// Defines HAS_UINT128 and uint128_t if applicable.
#include <stdbool.h>
#include "ryu/common.h"
#include "ryu/d2s_intrinsics.h"

static const uint64_t DOUBLE_POW5_INV_SPLIT2[15][2] = {
  {                    1u, 2305843009213693952u },
  {  5955668970331000884u, 1784059615882449851u },
  {  8982663654677661702u, 1380349269358112757u },
  {  7286864317269821294u, 2135987035920910082u },
  {  7005857020398200553u, 1652639921975621497u },
  { 17965325103354776697u, 1278668206209430417u },
  {  8928596168509315048u, 1978643211784836272u },
  { 10075671573058298858u, 1530901034580419511u },
  {   597001226353042382u, 1184477304306571148u },
  {  1527430471115325346u, 1832889850782397517u },
  { 12533209867169019542u, 1418129833677084982u },
  {  5577825024675947042u, 2194449627517475473u },
  { 11006974540203867551u, 1697873161311732311u },
  { 10313493231639821582u, 1313665730009899186u },
  { 12701016819766672773u, 2032799256770390445u }
};
static const uint32_t POW5_INV_OFFSETS[22] = {
  0x54544554, 0x04055545, 0x10041000, 0x00400414, 0x40010000, 0x41155555,
  0x00000454, 0x00010044, 0x40000000, 0x44000041, 0x50454450, 0x55550054,
  0x51655554, 0x40004000, 0x01000001, 0x00010500, 0x51515411, 0x05555554,
  0x50411500, 0x40040000, 0x05040110, 0x00000000,
};

static const uint64_t DOUBLE_POW5_SPLIT2[13][2] = {
  {                    0u, 1152921504606846976u },
  {                    0u, 1490116119384765625u },
  {  1032610780636961552u, 1925929944387235853u },
  {  7910200175544436838u, 1244603055572228341u },
  { 16941905809032713930u, 1608611746708759036u },
  { 13024893955298202172u, 2079081953128979843u },
  {  6607496772837067824u, 1343575221513417750u },
  { 17332926989895652603u, 1736530273035216783u },
  { 13037379183483547984u, 2244412773384604712u },
  {  1605989338741628675u, 1450417759929778918u },
  {  9630225068416591280u, 1874621017369538693u },
  {   665883850346957067u, 1211445438634777304u },
  { 14931890668723713708u, 1565756531257009982u }
};
static const uint32_t POW5_OFFSETS[21] = {
  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x40000000, 0x59695995,
  0x55545555, 0x56555515, 0x41150504, 0x40555410, 0x44555145, 0x44504540,
  0x45555550, 0x40004000, 0x96440440, 0x55565565, 0x54454045, 0x40154151,
  0x55559155, 0x51405555, 0x00000105
};

#define POW5_TABLE_SIZE 26
static const uint64_t DOUBLE_POW5_TABLE[POW5_TABLE_SIZE] = {
1ull, 5ull, 25ull, 125ull, 625ull, 3125ull, 15625ull, 78125ull, 390625ull,
1953125ull, 9765625ull, 48828125ull, 244140625ull, 1220703125ull, 6103515625ull,
30517578125ull, 152587890625ull, 762939453125ull, 3814697265625ull,
19073486328125ull, 95367431640625ull, 476837158203125ull,
2384185791015625ull, 11920928955078125ull, 59604644775390625ull,
298023223876953125ull //, 1490116119384765625ull
};

uint32_t __pow5Factor(uint64_t value) {
  const uint64_t m_inv_5 = 14757395258967641293u; // 5 * m_inv_5 = 1 (mod 2^64)
  const uint64_t n_div_5 = 3689348814741910323u;  // #{ n | n = 0 (mod 2^64) } = 2^64 / 5
  uint32_t count = 0;
  for (;;) {
    assert(value != 0);
    value *= m_inv_5;
    if (value > n_div_5)
      break;
    ++count;
  }
  return count;
}

#if defined(HAS_UINT128)

// Computes 5^i in the form required by Ryu, and stores it in the given pointer.
void __double_computePow5(const uint32_t i, uint64_t* const result) {
  const uint32_t base = i / POW5_TABLE_SIZE;
  const uint32_t base2 = base * POW5_TABLE_SIZE;
  const uint32_t offset = i - base2;
  const uint64_t* const mul = DOUBLE_POW5_SPLIT2[base];
  if (offset == 0) {
    result[0] = mul[0];
    result[1] = mul[1];
    return;
  }
  const uint64_t m = DOUBLE_POW5_TABLE[offset];
  const uint128_t b0 = ((uint128_t) m) * mul[0];
  const uint128_t b2 = ((uint128_t) m) * mul[1];
  const uint32_t delta = pow5bits(i) - pow5bits(base2);
  const uint128_t shiftedSum = (b0 >> delta) + (b2 << (64 - delta)) + ((POW5_OFFSETS[i / 16] >> ((i % 16) << 1)) & 3);
  result[0] = (uint64_t) shiftedSum;
  result[1] = (uint64_t) (shiftedSum >> 64);
}

// Computes 5^-i in the form required by Ryu, and stores it in the given pointer.
void __double_computeInvPow5(const uint32_t i, uint64_t* const result) {
  const uint32_t base = (i + POW5_TABLE_SIZE - 1) / POW5_TABLE_SIZE;
  const uint32_t base2 = base * POW5_TABLE_SIZE;
  const uint32_t offset = base2 - i;
  const uint64_t* const mul = DOUBLE_POW5_INV_SPLIT2[base]; // 1/5^base2
  if (offset == 0) {
    result[0] = mul[0];
    result[1] = mul[1];
    return;
  }
  const uint64_t m = DOUBLE_POW5_TABLE[offset]; // 5^offset
  const uint128_t b0 = ((uint128_t) m) * (mul[0] - 1);
  const uint128_t b2 = ((uint128_t) m) * mul[1]; // 1/5^base2 * 5^offset = 1/5^(base2-offset) = 1/5^i
  const uint32_t delta = pow5bits(base2) - pow5bits(i);
  const uint128_t shiftedSum =
    ((b0 >> delta) + (b2 << (64 - delta))) + 1 + ((POW5_INV_OFFSETS[i / 16] >> ((i % 16) << 1)) & 3);
  result[0] = (uint64_t) shiftedSum;
  result[1] = (uint64_t) (shiftedSum >> 64);
}

#else // defined(HAS_UINT128)

// Computes 5^i in the form required by Ryu, and stores it in the given pointer.
void __double_computePow5(const uint32_t i, uint64_t* const result) {
  const uint32_t base = i / POW5_TABLE_SIZE;
  const uint32_t base2 = base * POW5_TABLE_SIZE;
  const uint32_t offset = i - base2;
  const uint64_t* const mul = DOUBLE_POW5_SPLIT2[base];
  if (offset == 0) {
    result[0] = mul[0];
    result[1] = mul[1];
    return;
  }
  const uint64_t m = DOUBLE_POW5_TABLE[offset];
  uint64_t high1;
  const uint64_t low1 = umul128(m, mul[1], &high1);
  uint64_t high0;
  const uint64_t low0 = umul128(m, mul[0], &high0);
  const uint64_t sum = high0 + low1;
  if (sum < high0) {
    ++high1; // overflow into high1
  }
  // high1 | sum | low0
  const uint32_t delta = pow5bits(i) - pow5bits(base2);
  result[0] = shiftright128(low0, sum, delta) + ((POW5_OFFSETS[i / 16] >> ((i % 16) << 1)) & 3);
  result[1] = shiftright128(sum, high1, delta);
}

// Computes 5^-i in the form required by Ryu, and stores it in the given pointer.
void __double_computeInvPow5(const uint32_t i, uint64_t* const result) {
  const uint32_t base = (i + POW5_TABLE_SIZE - 1) / POW5_TABLE_SIZE;
  const uint32_t base2 = base * POW5_TABLE_SIZE;
  const uint32_t offset = base2 - i;
  const uint64_t* const mul = DOUBLE_POW5_INV_SPLIT2[base]; // 1/5^base2
  if (offset == 0) {
    result[0] = mul[0];
    result[1] = mul[1];
    return;
  }
  const uint64_t m = DOUBLE_POW5_TABLE[offset];
  uint64_t high1;
  const uint64_t low1 = umul128(m, mul[1], &high1);
  uint64_t high0;
  const uint64_t low0 = umul128(m, mul[0] - 1, &high0);
  const uint64_t sum = high0 + low1;
  if (sum < high0) {
    ++high1; // overflow into high1
  }
  // high1 | sum | low0
  const uint32_t delta = pow5bits(base2) - pow5bits(i);
  result[0] = shiftright128(low0, sum, delta) + 1 + ((POW5_INV_OFFSETS[i / 16] >> ((i % 16) << 1)) & 3);
  result[1] = shiftright128(sum, high1, delta);
}

#endif // defined(HAS_UINT128)


